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We investigate the steady-state ■probability density distribution of a 
large class of random processes by solving the governing Fokker-Planck 
equation. The random response statistics of a nonlinear single-degree-of- 
freedom mechanical model with hyperbolic tangent stiffness are discussed 
in some detail. The probability density of such systems is of the sech-power 
type which belongs to a class of distributions whose behaviors are carefully 
examined at the limits where the system parameter b approaches zero 
and infinity. Other important response statistics such as the mean square 
response, zero crossings, and peak distributions are also studied. 

I. INTRODUCTION 

In recent years, random vibrations of nonlinear systems have attracted 
considerable attention among engineers. 1 In this paper we investigate 
the Fokker-Planck equations 2,3 associated with a class of random proc- 
esses whose steady-state probability density distributions, of the Lia- 
punov potential function type. 

The random response statistics of a nonlinear single-degree-of -freedom 
model having a hyperbolic tangent stiffness function can be described 
as a softening spring whose force-deflection relationship is asymptotic 
to some maximum force level. Such a model can be used to represent 
an elastic-perfect-plastic system, material often encountered in classical 
mechanics. Limiting situations for a class of probability density func- 
tions such as those obtained in this study are examined. We show that 
the limiting behavior of the steady-state output probability density 
function of a system having a generalized hyperbolic tangent stiffness 
function, F(u) = (ko/b"' 1 ) tanh bu, is closely related to the range of the 
parameter a. At the limit b — > a> , the probability density function be- 
comes a Dirac delta (impulse) function or an exponential distribution, 
or identically approaches zero for all u, depending upon whether a 

543 



544 THE BELL SYSTEM TECHNICAL JOURNAL, APRIL 1970 

is less than, equal to, or greater than 1. At the limit b — * 0, it vanishes 
identically for all u and becomes a normal distribution or a Dirac delta 
function, depending upon whether a is less than, equal to, or greater 
than 2. In addition, we study statistics of other response parameters 
such as the mean square output, zero crossings and peak output dis- 
tribution, which are relevant to the control of the failure modes of the 
system. 

The motion of a dynamic system under purely random disturbance 
is described by a Markoff process y(0 = [yi(t) t y 2 (t), • • • , y n (t)] in the 
n-dimensional phase space. It can be shown 3,4 that for the initial con- 
dition 

P<J.) = II 5(2/. - Vi.) 

where y c is the initial state of y(0 and 8 is the Dirac delta function, the 
conditional probability density function p(y | y , t) of the process y(t) 
satisfies the forward Fokker-Planck equation, 

where 

GM = A,(y)p - * £ ■£■ [B u (j)p] (2) 

is the component of the probability current vector p(y\y,,t) in which 
AM = lim ( Vi ,, t - y<) (3) 

Al-«0 

and 

Ba(y) = lim ((y iiAl - y^iy,-,^ - y,)) (4) 

are intensity coefficients depending on the input and the properties of 
the system (the bracket indicating ensemble averaging). 

We are interested in the solution of the steady-state equation (1), 
that is when dp/dt = 0, for cases where all generalized response variables 
of a system in the 2?? phase-space coordinates are independent of one 
another. For this type of motion it is sometimes possible to find appro- 
priate partial operators which, when linearly operated on functions of the 
type Qiiydp + hi(^i)(dp/dy t ), generate an equation equivalent to (1). 
More specifically, the steady-state equation (1) can be put in the form 

ZL^gMp + hM^ =0 (5) 
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where the coefficients L, are arbitrary first-order partial-differential 
operators. If there exists a p(y) independent of initial conditions and 
satisfying each 

g<(y<)p + h t ( yi )&L = , (i= 1,2, ••• ,2n), 

then by Gray's uniqueness theorem such p(y) is the unique solution of 
equation (5). 8 Such a solution is 



P.i(y) = C II exp 



o hi(\i) 



(6) 



and C is the normalization factor. 

Equations (5) and (6) will be used in the following sections to analyze 
a class of nonlinear systems. 

II. HYPERBOLIC TANGENT STIFFNESS MODEL 

The mechanical system considered in this investigation is a single- 
degree-of-freedom oscillator with a mass m, a linear viscous damping 
c, and a nonlinear spring function F(u). When the system is subjected 
to a base acceleration excitation x b (t), its response is characterized by 
the displacement u(t) relative to the base. The equation of motion of 
the system is 

ii + 20u + ff(tt) = a(t) (7) 



where 



and 



2m m. 



a(t) = -x b (t). 



Let a(l) be a gaussian, stationary white noise with zero mean; that is, 
with the properties 

(a(t)) = 

(a(/,)a(/,)) = 2S„ S(U - t 2 ) 

where S is the constant power spectral density of a(t). Then the as- 
sociated steady-state Fokker-Planck equation for u(0 = [u(t), ii(t)} is 

^~ 2 p(n, u) - £ [up(u, u)] + ~ {[2/3* + S(u)]}p(u , u) = 0. (8) 
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For this two-dimensional case (n = 2), according to equations (5) 
and (6), the solution can be written down readily, 



p(«,*)-Cexp{-^ [J + />©#]} 

where C is the normalization factor determined by 



(9) 



// 



p(u, u) du dii = 1 . 



A special kind of softening spring described by a hyperbolic tangent 
function will now be considered. The force-deflection characteristic is 
shown in Fig. 1 and given as follows: 



k 
F(u) = -=r tanh bu, k , b > 0, 



(10) 



where k a is the initial stiffness, and b is the rate of convergence of the 
force-deflection curve. 

It should be noted that the spring force F(u) developed during the 
motion is bounded between kjb and — k /b. Therefore k /b may be 
regarded as yielding force and 1/6 the corresponding yielding displace- 
ment. The stiffness function F(u) described in equation (10) then pro- 
vides a good representation of the elastic-perfect-plastic behavior often 
encountered in the fields of classical mechanics and structural en- 
gineering. 

Let oil = kjm where u„ represents the natural frequency of the linear 
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Fig. 1 — Force-deflection relationship of hyperbolic tangent stiffness model. 
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oscillator with stiffness k ; then substitution of equation (10) into (9) 
yields 

p(u, ii) = C exp s — 2 2 — T2-2 hi cosh bu { (11) 

^ £0) o (T o <J ) 

where a 2 = (irS o /20ui 2 o ) is the variance of the linear response [that is, 
if F(u) = k a u]. 

Equation (11) shows that u and u are statistically independent. The 
probability density function for velocity u is normal with zero mean 
and variance <j 2 u> 2 , that is, 

P(u) = 7^-ti exp ( -5-3-2) , - °° < u < 00 . (12) 

The probability density function for the displacement u is 

p(u) = C 1 (6)[sech bu] 1Uo "" (13) 

where 

&(&) - (f sech 1/ffo36a btdtf \ (14) 

Because (sech bi-) i/b '' T ''' converges to zero very rapidly as £ — > 00 , 
Ci(b) in equation (14) can be evaluated numerically for any positive 
b. If l/b 2 cr 2 is an integer, equation (14) then becomes 

Cl(6) " 2-^7! B &D ~ 2k ~ X) (15) 

where 2D = l/b 2 a 2 „ are integers. 6 It is interesting to see that, if tanh bu 
is expanded into a power series, equation (13) then becomes 



p(u) = Cx(6) exp [-^2 (u 2 - |u 4 + • • •)] , 



1 ^ v 
U| =26' 



which indicates that a cubic softening spring with nonlinear coefficient 
k„b 2 /3 is the first approximation of the hyperbolic tangent spring. 
Values of p(u) given by equation (13) for various l/b 2 a 2 are shown in 
Figs. 2 and 3. 

III. LIMITING SITUATIONS OF p(u) 

In connection with the examination of the limiting behaviors of 
p(u) in equation (13), where the parameter b approaches zero and in- 
finity alternately, three useful theorems are presented. 
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-3 -2 



Fig. 2 — Sech-power probability density distributions. 




-s -6 



Fig. 3 — Sech-power probability density distributions. 
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Theorem 1: Let f n (x) be a sequence of nonnegative density functions 
integrable on [— °o , »]. Suppose there exists a sequence of positive inte- 
grable g n (x) such that 

g n (x) ^ F n (x) = /»(*)/£ /.(«) ds 
and 

lim / g n (x) dx + / g n (x) dx = for every c > 0. 

Then 

\imF„(x) = 8(x), (16) 

the Dirac delta function. 

Proof: We must show that, for every h in C°° (R), the space of test 
functions 

lim r F n (x)h(x)dx = h(0). 

n • r J — 00 

By the mean value theorem the following relationship holds: 
lim / F„(x)h(x) dx 

- lim f F n (x)h(x) dx + lim / F n (x)h(x) dx + lim h(£) 

n-»oo •» i n-«eo *'—oo n-«oo 

■/ F n {x)dx 

where ^ is some member of [— e, e], depending on e and n. The first two 
limits on the right side of the previous equation are zero by a comparison 
test; therefore, one can show that 

lim ( F n (x) dx = 1. 

n _,00 J-, 

Then 

lim hit) = lim j F n (x)h(x) dx. 

n—oo n—co J -no 

But the right side is independent of e. Thus, letting e approach zero, 



550 THE BELL SYSTEM TECHNICAL JOURNAL, APRIL 1970 

we deduce that 

ft(0) = lim f F n (x)h(x) dx, 

n-»oo •'-00 

which completes the proof. 

Instead of considering p(u) of equation (13), we shall investigate 
its more general form as 

[sech (bu)] 1/b " A 



P(6 .«>(«) = 



f [sech (by)] 1/b " A dy 

•'-00 



which is the steady-state displacement density function corresponding 
to a generalized hyperbolic tangent stiffness function 

F(u) = -%- tanh bu (17) 

if A = a\ . 

Theorem 2: Let p a (u) = linu^ Vtb. a )(u), then 
(i) a > 1 implies p a {u) = 0, 

(it) a = 1 implies p a (u) = (o~j) e ~ iuWA > 

and 

{Hi) a < 1 implies p a (u) = 5(w). 
Proof: First suppose a > 1. We observe that 

|P(».«)(«) 1 ^ ~p fora11 w > 

/ [sech (by)] 1 '*'* dy 

J — oo 

but 

[sech (by)] 1 ""* ^ exp (-| y |/6 a " l A), 
thus 

f 00 [sech (6z/)] 1/6ayl dy ^ f° exp (- 1 y \/b a ~ l A) dy = 2b'~ l A. 

J-09 "—08 

Thus, since 

I P( ». a) (u) | ^ l/ca-U) (is) 

for all u, we conclude that lima^ P(b, a )(u) = 0. 
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Now suppose that a = 1. Then 

lim In [sech (bu)] UbA = lim (1/bA) In (sech bu) 

6-»oo 6-«oo 

= lim (—1/bA) In (cosh bu) 

6->ao 

= lim [-tanh (bu)/A] = —\u\/A. 

o-»oo 

Thus, 

lim [sech (bu)] 1/bA = exp (- | u \/A). 

b—oa 

By the Lebesgue dominated-convergence theorem 

lim r [sech (by)] 1/bA dy = f exp (- 1 y \/A) dy = 2A. 

Thus, 

limp (6il) (w) = (2A)' 1 exp (-\u \/A). 

6-»oo 

Finally, we suppose that a < 1. Let 

2 1/i,aA exp(-|w I b a ~ a) /A) 
9a>.«M = — ! — ' 

j exp[-(6°- o) \y\)/A]dy 



or, equivalently, 

2A 



9<b. a) (u) = *— ^ exp (- | u | 6 (1 " a 7A). 



Then, since 

lim / g ib , a) (y) dy = lim / g ib . a) (y) dy = 

6-»oo ■» e i-»oo •'—oo 

for every e > and 

0c6.«)(w) ^ P(6.«>(w) for every w, 
we conclude from Theorem 1 that 

limp (6 , a) (u) = 5(w), 

b—co 

and the proof is completed. 

Theorem 3: Let p„(u) = lim ^ p<b. a )(u), ^en 
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(i) a > 2 implies p^{u) = 5(u), 

(it) a = 2 implies p' a (u) = normal distribution with variance a\ , and 

(Hi) a < 2 implies p«(u) = for all u. 

Proof: The case a > 2 implies pL(u) = 8{u) is proved in Appendix A 
in which we also show that there exists for every y e (0, 1) a C y > such 

that 



g b (u) = C y exp (_-^-)(26- I il) 



^_^wr_ for M <^. 

— /.«> ' ' 2d 

/ [sech MY' dy 

J -00 

It follows from the above that 

p a '(u) =0 for 1 < a < 2. 

From equation (18) we immediately have p'(u) = for a < 1. Now we 
have only to consider the cases when a = 2 and a = 1. In Appendix B 
we show that when a = 2, p„(«) is a normal distribution with variance 
o^ . In Appendix C we show that when a = 1, p„(w) = for all u. 
Therefore the proof is completed. 

According to Theorem 3, for p(u) given by equations (13) and (14), 
it follows that 

limp(w) = ,_ >i exp -~-g , a normal distribution, 

and according to Theorem 2 

limp(tt) = 0. (20) 

6-<oo 

At the limit 6 -* » , the yielding force fc /6 — » 0, that is, the system 
becomes perfect plastic. Thus one may expect an equal probability for 
all u on [— oo , oo ], as equation (20) indicated. As 6 — > 0, then kjb — > oo , 
the system remains elastic on [— oo , oo] with the initial stiffness k . It 
is well known that for linear systems the response probability distribu- 
tion is gaussian, which agrees with the result of equation (19). 

It is of interest to note that a similar force-deflection relationship as 
shown in Fig. 1 and as described by the hyperbolic tangent stiffness 
function given in equation (10) can be described by a full-wave smooth 
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limiter which is 



0(iO = h } Q exp (-Jp) dv 



in which d 2 = 2/irb 2 . 

It will be noted that in the above equation, G(u) is proportional to 
the integral of a gaussian probability curve. Function G(u) can also be 
used to evaluate the probability density function if made equivalent 
to F (u) as given by equation (10) when both G(u) and F(u) have the 
same initial slope and spring resistance limits. 

IV. OTHER IMPORTANT RESPONSE STATISTICS 

The failure modes of a mechanical system are generally controlled 
by response parameters such as the mean square displacement, zero 
crossings, or the peak displacement distributions. These response sta- 
tistics are closely related to p(u) and will be briefly discussed. 

The mean value of displacement response u vanishes because p(u) 
in equation (13) is an even function. The mean square or the variance 
of the displacement is given by 

a-Kjj) = (u 2 ) = J" u 2 p(u) du 

= 2d(6) [' u 2 sech 1/6 " ' bu du, (21) 

which can be evaluated in the following manner*: Let 

e ai dx 



J(a) = f 

J —a 



(cosh x) 2 " 
then it can be shown that 

and 

f ; x \ d \» = \\&> J(a) l = MWfW. 

J-oo (cosh x) L4 da J a=0 

where \p'(v) = (d/dv)[T'(v)/T(y)\ is the "trigamma" function and has 
been numerically tabulated. 7 



* This is pointed out by S. O. Rice. 
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From the above results and setting x = bu and v = bo in equation 
(21) we finally obtain 



o-l(jb) = ^ 3 C l (b)J(0)r(ba o ). 
Again according to Theorem 3, it is noted that 



(22) 



and from Theorem 2 that 



lim a 2 u (b) = o\ , 

b—0 



lim <rt(b) = co 



Thus the mean square response <r 2 u (b) with such limiting behavior can 
be illustrated as in Fig. 4. 

The expected number of zero crossings v + with positive slope per unit 
time (that is, the expected frequency) can be evaluated according to 
Rice, 8 



**00 

>* (b) = I up(0, u) du = 



Crib) (t 



2 \/3 



where Ci(b) is given by equation (14). 

Also according to Theorem 3, it can be shown that 

lim, + (6) =fS 

which is the frequency of the linear system. 



(23) 



(24) 



i/b 



Fig. 4 — Variation of mean-square displacement response. 
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The probability density of the peak amplitude of u(t), from equation 
(13), is given by 

dp(u) 



p(a) = - 



du 



[p(w)] u=0 



= 4r (sech ba) 1/a '"" tanh 6a. (25) 

(t o 

By the same argument used in the proof of Theorem 2 it can be 
shown that 

lim p(a) = 5(a), (26) 

6— » 

and it follows from Theorem 3 that 

. . ,. /tanh ba \ .. ,i/,.'t> i„ 

hmp(a) = lmi I tt-) lun sech ba 

6-0 6-0 \ 0- o 6 / t _ 



a I —a' 

= -2 exp 



(27) 



which is the Rayleigh distribution as expected because at this limit 
(b -> 0) the system becomes linear. The peak probability density dis- 
tribution p (a) for various b in equation (25) is illustrated in Fig. 5. 
Notice that for all cases p(a) approaches zero at large a; however, the 
rate of fall of p(a) is reduced as b is increased. 

It should be noted that when a = 1, the forcing function described 
by equation (17) approaches a sgn function as b — > », that is 

lim (k tanh bu) = k sgn u. 

Therefore, by taking appropriate limits to the density function pre- 
viously obtained for second-order systems with a general hyperbolic 
tangent forcing function, we obtain the steady-state solution for the 
response density of systems governed by the following equation 

u + 2/3u + fc sgn u = a(t). 

The response density for the above equation is given precisely in state- 
ment (ii) of Theorem 2, which can also be verified by using equation (9). 

V. ACKNOWLEDGMENTS 

The authors wish to thank S. O. Rice and M. Lax for their interest 
in this work and their many valuable comments and suggestions. 



556 THE BELL SYSTEM TECHNICAL JOURNAL, APRIL 1970 



0.5 


/ 




0.4 




\ b<T = 0, LINEAR CASE 


o 

b 
\ 

5. °- 3 




\ RALEIGH DISTRIBUTION 


Q. 






'* 






0.2 






0.1 






\ bo- = io 

1 1 1 1 1 — f. 



0.5 1.0 1.5 2.0 2.5 3.0 3.5 4.0 

a/<r 
Fig. 5 — Peak probability density distributions. 
APPENDIX A 

Partial Proof of Theorem 3 for Case a > 2 

We claim that if a > 2, then 

.. sech (bu) 1/b " A 

lim — - = 8(n). 



*° f sech (by) Ub " A dy 

•/ — on 



Proof: In view of Theorem 2 we have only to find functions 

sech (bu) 1/b " A 

a »( u ) ^ ~jz 

/ sech (by) 1/b " A dy 

J -00 

such that for every e > 

lim / g b (u) du = lim / g b (u) du = 

6-»0 «« fc-»0 «'-oo 

and such that sup /"^ g b (u) du < oo . We write 
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sech (by) l/baA = exp { - In [cosh (by)]/b a A \ . 
We observe that if | by | < ir/2, then 

In [cosh (by)]/b a A = E (l/6"A)(l/fc!)[lim (£)" In cosh (fiy)jb k 

= E (*//&" Afc!)[lim (£)' ' tanh (&,)]&'. 

We now make use of the fact that | by \ < 1 implies 

tanh (by) = ± /3„(- l) n+1 2 2 "(2 2 " - 1) &££-. 

Thus, there is for every y t (0, 1] a C\ > such that | bu \ < yir/2 implies 

sech (bu) 1/b ° A ^ C y exp (-< u 2 /b a ~ 2 A). 
Thus, since 

.. lexp (by) + exp (-&?/)] dy ~ 2 J \2 exp (by)) dy 

= 2b a ~ l A, 
we therefore take 

g b (u) = C y exp (-u 2 /b a ~ 2 A)2b a ~ 1 A for | u \ < yir/2b 
and 

g h (u) = 26"-^ sech (6M) 1/fca " 1 for \u\> yw/2b. 
Note that 

£ /11 «.W*.S26-Atg=)«4^P 

= 6—il E (fcvir) sech ™) 

< C 2 26 a_2 A. 
Thus, for a > 2 

lim P ff6 (u) dw = lim /" C 7 exp (-u 2 /b a - 2 A)2b a ~ 1 A 

X /i \l/b"A 

+ lim 26 0_2 A E (^Ttt) sech I™ 



dw 
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= lim^C 7 exp (-e 2 /b a - 2 A)2b a ~ l A + 

6-0 &0 



= 0. 

Also, we observe that 



lim / g b (u) du = lim / g b (u) du = 0. 

(,-,0 •'-oo b—Q ^ t 



APPENDIX B 



Partial Proof of Theorem 8 for Case a = 2 

Let p b (x) = 0(6)(sech bx) 1/b '"' , b ^ f or all x on [- °o , oo]. We will 
first show lim 6 - 0(b) = 0, then lim 6 ^ Pb(x) converges pointwisely to 
a normal distribution with zero mean and variance a\ . 

Proof: From the definition of Ptix), it follows that 

p b (x) = {exp [In 0(6)6Vj(sech bx)} 1/b "'°* 

i t \ i a(w i m sech bx 
\np b {x) = in 6{b) -\ =-5-3 

0" 

Then 

.. 1 . , , ,. —x tanh 6a; 

lim T2-2 In seen bx = lim kt~2 

6-0 o" 6-0 £ba 

by L'Hopital's rule. Thus, since 

,. tanh bx 

lim — =■ — = x, 

b-o O 

we conclude that 

x 2 
lim In p b (x) = —7T-2 , 

or that 

lim Pj, (x) = exp (— jr-g) , the normal distribution. 

By using these expressions and equation (17) we can conclude that 
Pa(u) is a normal distribution with variance <r\ when a = 2. 
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APPENDIX C 

Partial Proof of Theorem 3 for Case a = 1 
We claim that, when 

sech (bu) UbA 



P(6.1)(W) = ^ 

/ sech (by) UbA dy 

then 

li n P<6.i)('") = for all w. 

6-.0 

Proof: We observe that 

sech (6u) I/6il 



P(».i)(«) ^ 



f sech (6?/) 1/iJl <fy 

J_AT 



for all positive integers iV. 
We show that for every N > there is a ^ > such that b < n implies 

P(6.d(«) =5 1/22V. 

We can show, using L'Hopital's rule, that 

lim sech (bu) 1/bA = 1 

6-.0 

for all u and A. Thus, for every tj > no matter how small and every 
N' ^ N(l + r?)/(l — i]) we can show that there is a n > such that by 
taking b < n, we have 

P».i>(«) ^ * + ? 

/ (1 " V) dr, 
•'-If' 

Thus, b < n implies 

Thus, for every N > there is a /* > such that < b < /x implies 

P(6.d(«) ^ 2^- 
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Hence 

limp (6<1 ,(u) = 0. 

t-K) 

The proof is completed. 
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